/* Copyright (c) 2021, National University of Defense Technology. All rights reserved.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 * 
 *     http://www.apache.org/licenses/LICENSE-2.0
 * 
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include "parcsr_strength.h"

namespace YHAMG
{

#define ABS(x) ((x) > 0 ? (x) : -(x))

void ParCSR_Strength(const ParCSR_Matrix& A, bool* LocLoc_s, bool* LocExt_s, double threshold)
{
	int n = A.LocLoc.size[0];
	int* LocLoc_rowptr = A.LocLoc.rowptr;
	int* LocLoc_colind = A.LocLoc.colind;
	double* LocLoc_values = A.LocLoc.values;
	int* LocExt_rowptr = A.LocExt.rowptr;
	int* LocExt_colind = A.LocExt.colind;
	double* LocExt_values = A.LocExt.values;

#ifdef USE_OPENMP
#pragma omp parallel for schedule(guided)
#endif
	for (int i = 0; i < n; ++i)
	{
		double theta = 0.0;
		for (int j = LocLoc_rowptr[i]; j < LocLoc_rowptr[i + 1]; ++j)
		{
			if (LocLoc_colind[j] == i) continue;
			double abs_jval = ABS(LocLoc_values[j]);
			if (abs_jval > theta) theta = abs_jval;
		}
		for (int j = LocExt_rowptr[i]; j < LocExt_rowptr[i + 1]; ++j)
		{
			double abs_jval = ABS(LocExt_values[j]);
			if (abs_jval > theta) theta = abs_jval;
		}
		theta *= threshold;
		for (int j = LocLoc_rowptr[i]; j < LocLoc_rowptr[i + 1]; ++j)
			LocLoc_s[j] = LocLoc_colind[j] != i && ABS(LocLoc_values[j]) > theta;
		for (int j = LocExt_rowptr[i]; j < LocExt_rowptr[i + 1]; ++j)
			LocExt_s[j] = ABS(LocExt_values[j]) > theta;
	}
}

}